import os, itertools
import pandas as pd
import numpy as np
from bokeh.plotting import figure, output_file, show
from bokeh.io import output_notebook
from bokeh.palettes import Dark2_5, colorblind
from glob import glob
from bokeh.layouts import layout, column, row
from bokeh.models.widgets import Tabs, Panel
from bokeh.io import curdoc
from bokeh.plotting import figure
from scipy import interpolate, integrate
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
def read_magnitudes(fname):
COL_NAMES = ['time', 'U', 'B', 'V', 'R', 'I', 'J', 'H', 'K']
mag_table = pd.read_csv(fname, delim_whitespace=True, comment='#', names=COL_NAMES)
mag_table.iloc[:, 1:] = mag_table.iloc[:,1:].replace(99.999, np.nan)
return mag_table
def read_spectrum(fname):
col_names = ['wavelength']
with open(fname, encoding='utf8') as fh:
for i, line in enumerate(fh):
if i==2:
times = list(map(np.float64, line.strip().split(':')[1].split()))
spectrum = pd.read_csv(fname, delim_whitespace=True, comment='#', header=None, names=col_names + times)
return spectrum
def read_tgas(fname):
col_names = ['velocity']
with open(fname, encoding='utf8') as fh:
for i, line in enumerate(fh):
if i==2:
times = list(map(np.float64, line.strip().split(':')[1].split()))
spectrum = pd.read_csv(fname, delim_whitespace=True, comment='#', header=None, names=col_names + times)
return spectrum
# uncomment the following line to download the code comparison stuff
#!wget http://opensupernova.org/~wkerzend/files/ccompare1_models.tgz
#!tar zxf ccompare1_models.tgz
Run interactively on Google Colab on the following link
CODE_NAMES = ['artis', 'cmfgen', 'sedona', 'stella', 'sumo', 'supernu', 'urilight']
CODE_COLORS = {name:color for name, color in zip(CODE_NAMES, colorblind['Colorblind'][8])}
Integration of the synthetic spectra excluding gamma rays ($>50 Angstrom$)
output_notebook()
# create a new plot with a title and axis labels
fig_lbol = figure(title="Comparison of Lbol", x_axis_label='time [days]', y_axis_label='Log Lbol')
fig_edep = figure(title="Comparison of energy deposition", x_axis_label='time [days]', y_axis_label='Log Edeptot')
fig_econ = figure(title="Energy conservation", x_axis_label='time [days]', y_axis_label='(int t * lbol * dt)/(int t * edep * dt)')
lbol_data = {}
for code_name in CODE_NAMES:
if code_name == 'sedona':
fname = 'toy01/lbol_edep_toy01_{0}.txt'.format('sedonaexpansion')
else:
fname = 'toy01/lbol_edep_toy01_{0}.txt'.format(code_name)
if not os.path.exists(fname):
continue
lbol = pd.read_csv(fname, delim_whitespace=True, comment='#')
lbol.columns = ['time', 'lbol', 'edep']
lbol['log_lbol'] = np.log10(lbol.lbol)
lbol['log_edep'] = np.log10(lbol.edep)
lbol['econ']= np.hstack(([np.nan], integrate.cumtrapz(lbol.lbol*lbol.time, lbol.time) / integrate.cumtrapz(lbol.edep*lbol.time, lbol.time)))
fig_lbol.line(lbol.time, lbol.log_lbol, legend=code_name, line_width=2, color=CODE_COLORS[code_name])
fig_edep.line(lbol.time, lbol.log_edep, legend=code_name, line_width=2, color=CODE_COLORS[code_name])
fig_econ.line(lbol.time, lbol.econ, legend=code_name, line_width=2, color=CODE_COLORS[code_name])
lbol_data[code_name] = lbol
fig_lbol.legend.location = "top_right"
fig_lbol.legend.click_policy="hide"
fig_edep.legend.location = "top_right"
fig_edep.legend.click_policy="hide"
fig_econ.legend.location = "top_right"
fig_econ.legend.click_policy="hide"
# show the results
show(column([fig_lbol, fig_edep, fig_econ]))#
## Reading the data
tgas_data = {}
for code_name in CODE_NAMES:
if code_name == 'sedona':
fname = 'toy01/tgas_toy01_{0}.txt'.format('sedonaexpansion')
else:
fname = 'toy01/tgas_toy01_{0}.txt'.format(code_name)
if not os.path.exists(fname):
continue
try:
tgas_data[code_name] = read_tgas(fname)
except UnicodeDecodeError:
pass
except IndexError:
pass
TGAS_TIMES = [5, 10, 15, 20, 30, 50, 100, 200, 300, 400]
output_notebook()
interpolators = {}
for name in CODE_NAMES:
if name not in tgas_data:
continue
interpolators[name] = interpolate.interp1d(tgas_data[name].columns[1:],
tgas_data[name].iloc[:, 1:].values,
fill_value=np.nan, bounds_error=False)
all_tabs = []
for time in TGAS_TIMES:
fig = figure(x_axis_label='velocity [km/s]', y_axis_label='Tgas [K]')
for name, interpolator in interpolators.items():
velocity = tgas_data[name].velocity
fig.line(velocity, interpolator(time), legend=name, line_width=2, color=CODE_COLORS[name])
fig.legend.location = "top_right"
fig.legend.click_policy="hide"
lout = layout([[fig]])
tab = Panel(child=lout,title='t={0:.0f}'.format(time))
all_tabs.append(tab)
show(Tabs(tabs=all_tabs))
## Reading the data
mag_data = {}
for code_name in CODE_NAMES:
if code_name == 'sedona':
fname = 'toy01/mags_toy01_{0}.txt'.format('sedonaexpansion')
else:
fname = 'toy01/mags_toy01_{0}.txt'.format(code_name)
if not os.path.exists(fname):
continue
try:
mag_data[code_name] = read_magnitudes(fname)
except UnicodeDecodeError:
pass
except IndexError:
pass
output_notebook()
MAG_NAMES = list(mag_data.values())[0].columns[1:]
all_tabs = []
for mag in MAG_NAMES:
fig = figure()
for name, mag_table in mag_data.items():
fig.line(mag_table.time, mag_table[mag], legend=name, line_width=2, color=CODE_COLORS[name])
fig.y_range.flipped = True
fig.legend.location = "top_right"
fig.legend.click_policy="hide"
lout = layout([[fig]])
tab = Panel(child=lout,title=mag)
all_tabs.append(tab)
show(Tabs(tabs=all_tabs))
output_notebook()
MAG_NAMES = list(mag_data.values())[0].columns[1:]
all_tabs = []
for mag1, mag2 in zip(MAG_NAMES[:-1], MAG_NAMES[1:]):
fig = figure()
for name, mag_table in mag_data.items():
fig.line(mag_table.time, mag_table[mag1] - mag_table[mag2], legend=name, line_width=2, color=CODE_COLORS[name])
fig.y_range.flipped = True
fig.legend.location = "top_right"
fig.legend.click_policy="hide"
lout = layout([[fig]])
tab = Panel(child=lout,title='{0}-{1}'.format(mag1, mag2))
all_tabs.append(tab)
show(Tabs(tabs=all_tabs))
## Reading the data
spectral_data = {}
for code_name in CODE_NAMES:
if code_name == 'sedona':
fname = 'toy01/spectra_toy01_{0}.txt'.format('sedonaexpansion')
else:
fname = 'toy01/spectra_toy01_{0}.txt'.format(code_name)
if not os.path.exists(fname):
continue
try:
spectral_data[code_name] = read_spectrum(fname)
except UnicodeDecodeError:
pass
except IndexError:
pass
SPECTRAL_TIMES = [5, 10, 15, 20, 30, 50, 100, 200, 300, 400]
output_notebook()
interpolators = {}
for name in CODE_NAMES:
if name not in spectral_data:
continue
interpolators[name] = interpolate.interp1d(spectral_data[name].columns[1:],
spectral_data[name].iloc[:, 1:].values,
fill_value=np.nan, bounds_error=False)
all_tabs = []
for time in SPECTRAL_TIMES:
fig = figure()
for name, interpolator in interpolators.items():
wavelength = spectral_data[name].wavelength
fig.line(wavelength, interpolator(time), legend=name, line_width=2, color=CODE_COLORS[name])
fig.legend.location = "top_right"
fig.legend.click_policy="hide"
lout = layout([[fig]])
tab = Panel(child=lout,title='t={0:.0f}'.format(time))
all_tabs.append(tab)
show(Tabs(tabs=all_tabs))